# How Morality Politics determine Morality Policy Output ???
# Partisan Effects on Morality Policy Change

library(foreign)
library(interplot)
library(dummies)
library(logistf)
library(stargazer)
library(texreg)
library(dummies)
library(ggplot2)
library(ggpubr)

#Figure A1 (appendix): Number of reforms by policy area:

desc2 <- readRDS("data_reforms_by_policy.Rda")

jpeg("figure_A1_bw.jpg", width=4, height=4, units= 'in', res=600)
ggplot(desc2, aes(x=V1, y=number, fill=direction)) + 
  geom_bar(stat="identity", position=position_dodge(), color="black", width=1) +
  scale_fill_grey(start=0.8, end=0.2) +
  ylab("number of reforms") +
  xlab("policy area") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1.0))+
  scale_y_continuous(limits = c(0,50)) 
dev.off()

###########################################################################

#Figure A2 (appendix): Number of Reforms by year

ydata_all2 <- readRDS("data_reforms_by_year.Rda")

jpeg("figure_A2_bw.jpg", width=4, height=4, units= 'in', res=600)
ggplot(ydata_all2, aes(x=year, y=number, fill=direction )) + 
  geom_bar(stat="identity", position=position_dodge(), width=1) +
  scale_fill_manual(values = c("grey60", "grey0")) +
  ylab("number of reforms") +
  xlab("year") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1.0))+
  scale_y_continuous(limits = c(0,8), breaks = c(2,4,6,8) ) 
dev.off()

###############################################################################

#Figure 1 (main text): Number of Reforms by World

reli <- readRDS("data_agg_rel.Rda")
sec <- readRDS("data_agg_sec.Rda")

#construct individual parts of the graph (black and white and in colour)

p1 <- ggplot(reli, aes(x=country, y=number2, fill=direction)) + 
  geom_bar(stat="identity", position=position_dodge(), color="black", width=1) +
  scale_fill_grey() +
  ylab("number of reforms") +
  xlab("religious world")+
  theme(axis.text.x = element_text(angle = 90, vjust=0.5, hjust=1.0))+
  scale_y_continuous(limits = c(0,15), breaks=c(0,3,6,9,12,15)) 

p2 <-  ggplot(sec, aes(x=country, y=number2, fill=direction)) + 
  geom_bar(stat="identity", position=position_dodge(), color="black", width=1) +
  scale_fill_grey() +
  ylab("") +
  xlab("secular world")+
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1.0)) +
  scale_y_continuous(limits = c(0,15), breaks=c(0,3,6,9,12,15))

# combine individual graphs:  

jpeg("figure_1_bw.jpg", width=4, height=4, units='in', res=600)
ggarrange(p1, p2, 
          ncol = 2, nrow=1,align="h", common.legend = TRUE, legend = "bottom")
dev.off()  

#################################################################################
# customized extract function for texreg
extract.logistf <- function(model){
  s <- summary(model)
  names <- s$terms
  co <- s$coefficients
  se <- sqrt(diag(s$var))
  pval <- s$prob
  n <- nobs(model)
  eaic <- extractAIC(s)
  eaic2 <- eaic[2]
  gof <- c(eaic2, n)
  gof.names <- c("AIC","Obs.")
  
  tr <- createTexreg(  
    coef.names = names, 
    coef = co,
    se = se,
    pvalues = pval,
    gof.names = gof.names,
    gof = gof
  )
  return(tr)  
}    

setMethod("extract", signature = className("logistf", "logistf"), definition = extract.logistf)

##################################################################################
# Regression Analysis
      policy_sec <- readRDS("policy_sec.Rda")

      # Secular World: Liberalization (H1a) TABLE A1:
                #Table A1 - model 1:
                    sec_lib_1c <- logistf(data=policy_sec, dvlibreform ~ cons_ca + v7210_pii + kat + loggdb + tsR2 + i_all_perc 
                      , firth=TRUE)
                    
                    summary(sec_lib_1c)

                #Table A1 - model 2:  
                    sec_lib_2c <- logistf(data=policy_sec, dvlibreform ~ cons_ca + v7210_pii + kat + loggdb + tsR2 + i_all_perc 
                                          #   + policy_sec_1 # abortion
                                          + policy_sec_2 # euthanasia
                                          + policy_sec_3 # homo
                                          + policy_sec_4 # pornography
                                          + policy_sec_5 # prostitution
                                          + policy_sec_6 # ssm
                                          , firth=TRUE)
                    summary(sec_lib_2c)
                    
                #Table A1 - model 3 
                    sec_lib_3c <- logistf(data=policy_sec, dvlibreform ~ cons_ca + v7210_pii + kat + loggdb + tsR2 + i_all_perc 
                                          #   + policy_sec$policy_sec_1 # abortion
                                          + policy_sec_2 # euthanasia
                                          + policy_sec_3 # homo
                                          + policy_sec_4 # pornography
                                          + policy_sec_5 # prostitution
                                          + policy_sec_6 # ssm
                                          + tsR_abortion2
                                          + tsR_euth2
                                          + tsR_homo2
                                          + tsR_porn2
                                          + tsR_pro2
                                          + tsR_ssm2
                                          , firth=TRUE)
                    summary(sec_lib_3c)

                #Table A 1 - model 4
                    sec_lib_4c <- logistf(data=policy_sec, dvlibreform ~ cons_ca + v7210_pii + kat + loggdb + tsR2 + i_all_perc 
                                          #   + sixties
                                          + seventies
                                          + eighties
                                          + ninties
                                          + nuller
                                          , firth=TRUE)
                    summary(sec_lib_4c)
               
              # Table A 1 - model 5     
                    sec_lib_5c <- logistf(data=policy_sec, dvlibreform ~ cons_ca + v7210_pii + kat + loggdb + tsR2 + i_all_perc 
                                          + sec_6
                                          + sec_7 
                                          + sec_8 
                                          + sec_9 
                                          + sec_18
                                          + sec_21
                                          #    + sec_26
                                          , firth=TRUE)
                    summary(sec_lib_5c) 
                    #Secular World: Liberalization (H1):

                    texreg(list(sec_lib_1c, sec_lib_2c, sec_lib_3c, sec_lib_4c, sec_lib_5c),
                         #   file="tableA1.doc",
                            custom.coef.names = c(
                             "Intercept",
                             "Strength Conservative Parties",
                             "Political constraints",
                             "Catholicism",
                             "GDP/capita (log)",
                             "time since last reform",
                             "level of permissiveness",
                             "euthansia",
                             "homosexuality",
                             "pornography",
                             "prostitution",
                             "same-sex marriage",
                             "time since reform (abortion)",
                             "time since reform (euthanasia)",
                             "time since reform (homosexuality)",
                             "time since reform (pornography)",
                             "time since reform (prostitution)",
                             "time since reform (same-sex marriage)",
                             "70s",
                             "80s",
                             "90s",
                             "2000s",
                             "Denmark",
                             "Finland",
                             "France",
                             "Greece",
                             "Norway",
                             "Sweden"
                           ),
                           reorder.coef = c(2:28,1),
                           caption= "Secular World: Liberalization (H1)",
                         #  label="tab:2",
                           single.row=TRUE,
                           sideways=TRUE
                    )
 
                    
    #############################################################################                  
        #Secular World - Restrictive Reforms (H1b) - Table A2.1
              #Table 2.1:
                    #Model 1 - 5
                    sr1 <- logistf(data=policy_sec, dvresreform ~ cons_ca
                                   , firth=TRUE)
                    summary(sr1)
                    #
                    sr2 <- logistf(data=policy_sec, dvresreform ~ cons_ca + v7210_pii
                                   , firth=TRUE)
                    summary(sr2)
                    #
                    sr3 <- logistf(data=policy_sec, dvresreform ~ cons_ca + kat
                                   , firth=TRUE)
                    summary(sr3)
                    #
                    sr4 <- logistf(data=policy_sec, dvresreform ~ cons_ca + loggdb
                                   , firth=TRUE)
                    summary(sr4)
                    #
                    sr5 <- logistf(data=policy_sec, dvresreform ~ cons_ca + tsR2
                                   , firth=TRUE)
                    summary(sr5)
                    #
                    sr6 <- logistf(data=policy_sec, dvresreform ~ cons_ca + i_all_perc
                                   , firth=TRUE)
                    summary(sr6)
                    #
                    sr7 <- logistf(data=policy_sec, dvresreform ~ cons_ca + policy_sec_1
                                   , firth=TRUE)
                    summary(sr7)
                    #
                    sr8 <- logistf(data=policy_sec, dvresreform ~ cons_ca + policy_sec_2
                                   , firth=TRUE)
                    summary(sr8)
                    #
                    sr9 <- logistf(data=policy_sec, dvresreform ~ cons_ca + policy_sec_3
                                   , firth=TRUE)
                    summary(sr9)
                    #
                    sr10 <- logistf(data=policy_sec, dvresreform ~ cons_ca + policy_sec_4
                                    , firth=TRUE)
                    summary(sr10)
                    #
                    sr11 <- logistf(data=policy_sec, dvresreform ~ cons_ca + policy_sec_5
                                    , firth=TRUE)
                    summary(sr11)
                    #
                    sr12 <- logistf(data=policy_sec, dvresreform ~ cons_ca + policy_sec_6
                                    , firth=TRUE)
                    summary(sr12)
                    #
                    sr13 <- logistf(data=policy_sec, dvresreform ~ cons_ca + tsR_abortion2
                                    , firth=TRUE)
                    summary(sr13)
                    #
                    sr14 <- logistf(data=policy_sec, dvresreform ~ cons_ca + tsR_euth2
                                    , firth=TRUE)
                    summary(sr14)
                    #
                    sr15 <- logistf(data=policy_sec, dvresreform ~ cons_ca + tsR_homo2
                                    , firth=TRUE)
                    summary(sr15)
                    #
                    sr16 <- logistf(data=policy_sec, dvresreform ~ cons_ca + tsR_porn2
                                    , firth=TRUE)
                    summary(sr16)
                    #
                    sr17 <- logistf(data=policy_sec, dvresreform ~ cons_ca + tsR_pro2
                                    , firth=TRUE)
                    summary(sr17)
                    #
                    sr18 <- logistf(data=policy_sec, dvresreform ~ cons_ca + tsR_ssm2
                                    , firth=TRUE)
                    summary(sr18)
                    #
                    sr19 <- logistf(data=policy_sec, dvresreform ~ cons_ca + sixties
                                    , firth=TRUE)
                    summary(sr19)
                    #
                    sr20 <- logistf(data=policy_sec, dvresreform ~ cons_ca + seventies
                                    , firth=TRUE)
                    summary(sr20)
                    #
                    sr21 <- logistf(data=policy_sec, dvresreform ~ cons_ca + eighties
                                    , firth=TRUE)
                    summary(sr21)
                    #
                    sr22 <- logistf(data=policy_sec, dvresreform ~ cons_ca + ninties
                                    , firth=TRUE)
                    summary(sr22)
                    #
                    sr23 <- logistf(data=policy_sec, dvresreform ~ cons_ca + nuller
                                    , firth=TRUE)
                    summary(sr23)
                    #
                    sr24b <- logistf(data=policy_sec, dvresreform ~ cons_ca + sec_6
                                     , firth=TRUE)
                    summary(sr24b)
                    
                    #
                    sr24 <- logistf(data=policy_sec, dvresreform ~ cons_ca + sec_7
                                    , firth=TRUE)
                    summary(sr24)
                    #
                    sr25 <- logistf(data=policy_sec, dvresreform ~ cons_ca + sec_8
                                    , firth=TRUE)
                    summary(sr25)
                    #
                    sr26 <- logistf(data=policy_sec, dvresreform ~ cons_ca + sec_9
                                    , firth=TRUE)
                    summary(sr26)
                    #
                    sr27 <- logistf(data=policy_sec, dvresreform ~ cons_ca + sec_18
                                    , firth=TRUE)
                    summary(sr27)
                    #
                    sr28 <- logistf(data=policy_sec, dvresreform ~ cons_ca + sec_21
                                    , firth=TRUE)
                    summary(sr28)
                    #
                    sr29 <- logistf(data=policy_sec, dvresreform ~ cons_ca + sec_26
                                    , firth=TRUE)
                    summary(sr29)                    
                    
                    
                    
                    
                    
                    # Table A2.1:
                    texreg(list(sr1, sr2, sr3, sr4, sr5),
                       #    file="tableA2.1.doc",
                           custom.coef.names = c(
                             "Intercept",
                             "Strength Conservative Parties",
                             "Political constraints",
                             "Catholicism",
                             "GDP/capita (log)",
                             "time since last reform"
                           ),
                           reorder.coef = c(2:6,1),
                           caption= "Secular World: Restrictive Reform (H2)",
                    #       label="tab:3",
                           single.row=TRUE
                    #       sideways=TRUE
                    )
                    
                    
                    # Table A2.2:
                    texreg(list(sr6, sr7, sr8, sr9, sr10),
                        #   file="tableA2.2.doc",
                           custom.coef.names = c(
                             "Intercept",
                             "Strength Conservative Parties",
                             "level of permissiveness",
                             "abortion",
                             "euthansia",
                             "homosexuality",
                             "pornography"
                           ),
                           reorder.coef = c(2:7,1),
                           caption= "Secular World: Restrictive Reform (H2)",
                     #      label="tab:3",
                           single.row=TRUE
                    #       sideways=TRUE
                    )
                    
                    # Table A2.3:
                    texreg(list(sr11, sr12, sr13, sr14, sr15),
                        #   file="tableA2.3.doc",
                           custom.coef.names = c(
                             "Intercept",
                             "Strength Conservative Parties",
                             "prostitution",
                             "same-sex marriage",
                             "time since reform (abortion)",
                             "time since reform (euthanasia)",
                             "time since reform (homosexuality)"
                           ),
                           reorder.coef = c(2:7,1),
                           caption= "Secular World: Restrictive Reform (H2)",
                      #     label="tab:3",
                           single.row=TRUE
                      #     sideways=TRUE
                    )
                    
                    # Table A2.4:
                    texreg(list(sr16, sr17, sr18, sr19, sr20),
                         #  file="tableA2.4.doc",
                           custom.coef.names = c(
                             "Intercept",
                             "Strength Conservative Parties",
                             "time since reform (pornography)",
                             "time since reform (prostitution)",
                             "time since reform (same-sex marriage)",
                             "60s",
                             "70s"
                           ),
                           reorder.coef = c(2:7,1),
                           caption= "Secular World: Restrictive Reform (H2)",
                      #     label="tab:3",
                           single.row=TRUE
                      #     sideways=TRUE
                    )
                    
                    # Table A2.5:
                    texreg(list(sr21, sr22, sr23, sr24b, sr24),
                         #  file="tableA2.5.doc",
                           custom.coef.names = c(
                             "Intercept",
                             "Strength Conservative Parties",
                             "80s",
                             "90s",
                             "2000s",
                             "Denmark",
                             "Finland"
                          ),
                           reorder.coef = c(2:7,1),
                           caption= "Secular World: Restrictive Reform (H2)",
                      #     label="tab:3",
                           single.row=TRUE
                      #     sideways=TRUE
                    )
                    
                    
                    # Table A2.6:
                    texreg(list(sr25, sr26, sr27, sr28, sr29),
                         #  file="tableA2.6.doc",
                           custom.coef.names = c(
                             "Intercept",
                             "Strength Conservative Parties",
                             "France",
                             "Greece",
                             "Norway",
                             "Sweden",
                             "United Kingdom"
                           ),
                           reorder.coef = c(2:7,1),
                           caption= "Secular World: Restrictive Reform (H2)",
                     #      label="tab:3",
                           single.row=TRUE
                     #      sideways=TRUE
                    )
##################################################
          
            policy_rel <- readRDS("policy_rel.Rda")
            
                # Tables A3: Secular World Restrictive Reform: Models 1 -33:

                  rr1 <- logistf(data=policy_rel, dvresreform ~ cons_ca
                                   , firth=TRUE)
                    summary(rr1)
                    #
                    rr2 <- logistf(data=policy_rel, dvresreform ~ cons_ca + v7210_pii
                                   , firth=TRUE)
                    summary(rr2)
                    #
                    rr3 <- logistf(data=policy_rel, dvresreform ~ cons_ca + kat
                                   , firth=TRUE)
                    summary(rr3)
                    #
                    rr4 <- logistf(data=policy_rel, dvresreform ~ cons_ca + loggdb
                                   , firth=TRUE)
                    summary(rr4)
                    #
                    rr5 <- logistf(data=policy_rel, dvresreform ~ cons_ca + tsR2
                                   , firth=TRUE)
                    summary(rr5)
                    #
                    rr6 <- logistf(data=policy_rel, dvresreform ~ cons_ca + i_all_perc
                                   , firth=TRUE)
                    summary(rr6)
                    #
                    rr7 <- logistf(data=policy_rel, dvresreform ~ cons_ca + policy_rel_1
                                   , firth=TRUE)
                    summary(rr7)
                    #
                    rr8 <- logistf(data=policy_rel, dvresreform ~ cons_ca + policy_rel_2
                                   , firth=TRUE)
                    summary(rr8)
                    #
                    rr9 <- logistf(data=policy_rel, dvresreform ~ cons_ca + policy_rel_3
                                   , firth=TRUE)
                    summary(rr9)
                    #
                    rr10 <- logistf(data=policy_rel, dvresreform ~ cons_ca + policy_rel_4
                                    , firth=TRUE)
                    summary(rr10)
                    #
                    rr11 <- logistf(data=policy_rel, dvresreform ~ cons_ca + policy_rel_5
                                    , firth=TRUE)
                    summary(rr11)
                    #
                    rr12 <- logistf(data=policy_rel, dvresreform ~ cons_ca + policy_rel_6
                                    , firth=TRUE)
                    summary(rr12)
                    #
                    rr13 <- logistf(data=policy_rel, dvresreform ~ cons_ca + tsR_abortion2
                                    , firth=TRUE)
                    summary(rr13)
                    #
                    rr14 <- logistf(data=policy_rel, dvresreform ~ cons_ca + tsR_euth2
                                    , firth=TRUE)
                    summary(rr14)
                    #
                    rr15 <- logistf(data=policy_rel, dvresreform ~ cons_ca + tsR_homo2
                                    , firth=TRUE)
                    summary(rr15)
                    #
                    rr16 <- logistf(data=policy_rel, dvresreform ~ cons_ca + tsR_porn2
                                    , firth=TRUE)
                    summary(rr16)
                    #
                    rr17 <- logistf(data=policy_rel, dvresreform ~ cons_ca + tsR_pro2
                                    , firth=TRUE)
                    summary(rr17)
                    #
                    rr18 <- logistf(data=policy_rel, dvresreform ~ cons_ca + tsR_ssm2
                                    , firth=TRUE)
                    summary(rr18)
                    #
                    rr19 <- logistf(data=policy_rel, dvresreform ~ cons_ca + sixties
                                    , firth=TRUE)
                    summary(rr19)
                    #
                    rr20 <- logistf(data=policy_rel, dvresreform ~ cons_ca + seventies
                                    , firth=TRUE)
                    summary(rr20)
                    #
                    rr21 <- logistf(data=policy_rel, dvresreform ~ cons_ca + eighties
                                    , firth=TRUE)
                    summary(rr21)
                    #
                    rr22 <- logistf(data=policy_rel, dvresreform ~ cons_ca + ninties
                                    , firth=TRUE)
                    summary(rr22)
                    #
                    rr23 <- logistf(data=policy_rel, dvresreform ~ cons_ca + nuller
                                    , firth=TRUE)
                    summary(rr23)
                    #
                    rr24 <- logistf(data=policy_rel, dvresreform ~ cons_ca + rel_1
                                    , firth=TRUE)
                    summary(rr24)
                    #
                    rr25 <- logistf(data=policy_rel, dvresreform ~ cons_ca + rel_3
                                    , firth=TRUE)
                    summary(rr25)
                    #
                    #
                    rr27 <- logistf(data=policy_rel, dvresreform ~ cons_ca + rel_10
                                    , firth=TRUE)
                    summary(rr27)
                    #
                    rr28 <- logistf(data=policy_rel, dvresreform ~ cons_ca + rel_12
                                    , firth=TRUE)
                    summary(rr28)
                    #
                    rr29 <- logistf(data=policy_rel, dvresreform ~ cons_ca + rel_13
                                    , firth=TRUE)
                    summary(rr29)
                    
                    rr30 <- logistf(data=policy_rel, dvresreform ~ cons_ca + rel_17
                                    , firth=TRUE)
                    summary(rr30)
                    
                    rr31 <- logistf(data=policy_rel, dvresreform ~ cons_ca + rel_19
                                    , firth=TRUE)
                    summary(rr31)
                    
                    rr32 <- logistf(data=policy_rel, dvresreform ~ cons_ca + rel_22
                                    , firth=TRUE)
                    summary(rr32)
                    
                    rr33 <- logistf(data=policy_rel, dvresreform ~ cons_ca + rel_23
                                    , firth=TRUE)
                    summary(rr33)
                    
      ###################################
                    
                    texreg(list(rr1, rr2, rr3, rr4, rr5),
                         #   file="tableA3.1.doc",
                           custom.coef.names = c(
                             "Intercept",
                             "Strength Conservative Parties",
                             "Political constraints",
                             "Catholicism",
                             "GDP/capita (log)",
                             "time since last reform"
                           ),
                           reorder.coef = c(2:6,1),
                           caption= "Religious World: Restrictive Reform (H4)",
                    #       label="tab:5",
                           single.row=TRUE,
                     #      sideways=TRUE
                    )
                      ####################################
                    
                    texreg(list(rr6, rr7, rr8, rr9, rr10),
                       #    file="tableA3.2.doc",
                           custom.coef.names = c(
                             "Intercept",
                             "Strength Conservative Parties",
                             "level of permissiveness",
                             "abortion",
                             "euthanasia",
                             "homosexuality",
                             "pornography"
                           ),
                           reorder.coef = c(2:7,1),
                           caption= "Religious World: Restrictive Reform (H4)",
                      #     label="tab:5",
                           single.row=TRUE
                      #     sideways=TRUE
                    )
  
                    ############################
                    texreg(list(rr11, rr12, rr13, rr14, rr15),
                       #    file="tableA3.3.doc",
                           custom.coef.names = c(
                             "Intercept",
                             "Strength Conservative Parties",
                             "prostitution",
                             "same-sex marriage",
                             "time since reform (abortion)",
                             "time since reform (euthanasia)",
                             "time since reform (homosexuality)"
                           ),
                           reorder.coef = c(2:7,1),
                           caption= "Religious World: Restrictive Reform (H4)",
                     #      label="tab:5",
                           single.row=TRUE
                      #     sideways=TRUE
                    )
                    ############################
                    texreg(list(rr16, rr17, rr18, rr19, rr20),
                        #  file="tableA3.4.doc",
                               custom.coef.names = c(
                             "Intercept",
                             "Strength Conservative Parties",
                             "time since reform (pornography)",
                             "time since reform (prostitution)",
                             "time since reform (same-sex marriage)",
                             "60s",
                             "70s"
                           ),
                           reorder.coef = c(2:7,1),
                           caption= "Religious World: Restrictive Reform (H4)",
                   #        label="tab:5",
                           single.row=TRUE
                  #         sideways=TRUE
                    )
                    
                    ############################
                    texreg(list(rr21, rr22, rr23, rr24, rr25),
                         #   file="tableA3.5.doc",
                           custom.coef.names = c(
                             "Intercept",
                             "Strength Conservative Parties",
                             "80s",
                             "90s",
                             "2000s",
                             "Austria",
                             "Belgium"
                           ),
                           reorder.coef = c(2:7,1),
                           caption= "Religious World: Restrictive Reform (H4)",
                     #      label="tab:5",
                           single.row=TRUE
                     #      sideways=TRUE
                    )
                    
                    ############################
                    texreg(list(rr27, rr28, rr29, rr30),
                      #    file = "tableA3.6.doc",
                             custom.coef.names = c(
                             "Intercept",
                             "Strength Conservative Parties",
                             "Germany",
                             "Ireland",
                             "Italy",
                             "Netherlands"
                           ),
                           reorder.coef = c(2:6,1),
                           caption= "Religious World: Restrictive Reform (H4)",
                       #    label="tab:5",
                           single.row=TRUE
                       #    sideways=TRUE
                    )
                    
                    ##############################
                    texreg(list(rr31, rr32, rr33),
                         #   file="tableA3.7.doc",
                           custom.coef.names = c(
                             "Intercept",
                             "Strength Conservative Parties",
                             "Portugal",
                             "Spain",
                             "Switzerland"
                           ),
                           reorder.coef = c(2:5,1),
                           caption= "Religious World: Restrictive Reform (H4)",
                      #     label="tab:5",
                           single.row=TRUE
                      #     sideways=TRUE
                    )
                    
                    ###############################
    #H3: Religious World - Liberalization Table A4 models 1-5:      

                    rel_lib_1c <- logistf(data=policy_rel, dvlibreform ~ cons_ca + v7210_pii + kat + loggdb + tsR2 + i_all_perc 
                                          , firth=TRUE, method="Penalized ML")
                    summary(rel_lib_1c)

                    #
                    rel_lib_2c <- logistf(data=policy_rel, dvlibreform ~ cons_ca + v7210_pii + kat + loggdb + tsR2 + i_all_perc 
                                          #       + policy_rel_1      ##    reform is a abortion reform
                                          + policy_rel_2      ##    reform is a euthanasia reform  
                                          + policy_rel_3      ##    reform is a homosexuality reform
                                          + policy_rel_4      ##    reform is a pornography reform
                                          + policy_rel_5      ##    reform is a prostitution reform
                                          + policy_rel_6      ##    reform is a ssm reform
                                          , firth=TRUE)
                    summary(rel_lib_2c)
                    
                   #
                    rel_lib_3c <- logistf(data=policy_rel, dvlibreform ~ cons_ca + v7210_pii + kat + loggdb + tsR2 + i_all_perc 
                                          #       + policy_rel_1      ##    reform is a abortion reform
                                          + policy_rel_2      ##    reform is a euthanasia reform  
                                          + policy_rel_3      ##    reform is a homosexuality reform
                                          + policy_rel_4      ##    reform is a pornography reform
                                          + policy_rel_5      ##    reform is a prostitution reform
                                          + policy_rel_6      ##    reform is a ssm reform
                                          + tsR_abortion2
                                          + tsR_euth2
                                          + tsR_homo2
                                          + tsR_porn2
                                          + tsR_pro2
                                          + tsR_ssm2
                                          , firth=TRUE)
                    summary(rel_lib_3c)
                    #4-L
                   #
                    rel_lib_4c <- logistf(data=policy_rel, dvlibreform ~ cons_ca + v7210_pii + kat + loggdb + tsR2 + i_all_perc 
                                          #   + sixties
                                          + seventies
                                          + eighties
                                          + ninties
                                          + nuller
                                          , firth=TRUE)
                    summary(rel_lib_4c) 
                     #
                    rel_lib_5c <- logistf(data=policy_rel, dvlibreform ~ cons_ca + v7210_pii + kat + loggdb + tsR2 + i_all_perc 
                                          + rel_1 
                                          + rel_3
                                    #      + rel_6
                                          #   + rel_10 
                                          + rel_12
                                          + rel_13 
                                          + rel_17 
                                          + rel_19
                                          + rel_22
                                          + rel_23
                                          , firth=TRUE)
                    summary(rel_lib_5c) 
                    

                    
                    
                    
                    
  ############     
                    texreg(list(rel_lib_1c, rel_lib_2c, rel_lib_3c, rel_lib_4c, rel_lib_5c),
                        #   file="tableA4.doc",
                           custom.coef.names = c(
                             "Intercept",
                             "Strength Conservative Parties",
                             "Political constraints",
                             "Catholicism",
                             "GDP/capita (log)",
                             "time since last reform",
                             "level of permissiveness",
                             "euthanasia",
                             "homosexuality",
                             "pornography",
                             "prostitution",
                             "same-sex marriage",
                             "time since reform (abortion)",
                             "time since reform (euthanasia)",
                             "time since reform (homosexuality)",
                             "time since reform (pornography)",
                             "time since reform (prostitution)",
                             "time since reform (same-sex marriage)",
                             "70s",
                             "80s",
                             "90s",
                             "2000s",
                             "Austria",
                             "Belgium",
                        #     "Denmark",
                             "Ireland",
                             "Italy",
                             "Netherlands",
                             "Portugal",
                             "Spain",
                             "Switzerland"
                           ),
                           reorder.coef = c(2:30,1),
                           caption= "Religious World: Liberalization (H3)",
                      #     label="tab:4",
                           single.row=TRUE
                      #     sideways=TRUE
                    )

                    
####################################################    
 #main figure regression output (main text)
                    library(forestplot)
                    moral <- structure(list(
                      mean = c(sec_lib_1c$coefficients[2], sr1$coefficients[2], rr1$coefficients[2], rel_lib_1c$coefficients[2]),
                      lower = c(sec_lib_1c$ci.lower[2], sr1$ci.lower[2], rr1$ci.lower[2], rel_lib_1c$ci.lower[2]),
                      upper = c(sec_lib_1c$ci.upper[2], sr1$ci.upper[2], rr1$ci.upper[2], rel_lib_1c$ci.upper[2])),
                      .Names = c("mean", "lower", "upper"),
                      row.names = c(NA,-4L),
                      class ="data.frame")
                    
                    tabletext2 <- c("H1a: Secular World - Liberalization", "H1b: Secular World - Restrictive Reform", "H2: Religious World - Restrictive Reform", "H3: Religious World - Liberalization")
                    
                    jpeg("figure_2_bw.jpg", width=6, height=4, units='in', res=600)
                    
                    forestplot(tabletext2, moral, new_page= FALSE,
                               zero=0, lwd.zero=1,
                               txt_gp=fpTxtGp(label = list(gpar(fontsize = 12)), ticks=gpar(fontsize=12, cex=1),
                                              xlab=gpar(fontsize=12, cex=1)),
                               col=fpColors(lines="black", box="black", zero="black"),
                               boxsize=0.05, clip=c(-Inf,Inf), xticks=c(-5,-4,-3,-2,-1,0,1,2,3,4,5)
                    )
                    
                    dev.off()
                    

####################################################
  # Additional Robustness Checks:
  # exclude Ireland and Italy from religious world and add them to the secular world
  
  # Table A5.1: Religious World without Italy and Ireland:
  rr1_test <- logistf(data=subset(policy_rel,c_id < 12| c_id > 13), dvresreform ~ cons_ca
                                      , firth=TRUE)
  summary(rr1_test)
  #
  rel_lib_1c_test <- logistf(data=subset(policy_rel,c_id < 12| c_id > 13), dvlibreform ~ cons_ca
                                               , firth=TRUE)
  summary(rr1_test)
            
  texreg(list(rr1_test, rel_lib_1c_test),
     #     file="tableA5.1.doc",
         custom.coef.names = c(
           "Intercept",
           "Strength Conservative Parties"
         ),
         reorder.coef = c(2,1),
         caption= "Religious World: without Italy and Ireland",
      #   label="tab:5",
         single.row=TRUE
      #   sideways=TRUE
  )
  
  # Table A5.2: Secular World with Italy and Ireland:
  policy_sec_robust <- readRDS("secITAIRE.Rda")
  #
  sr1_test <- logistf(data=policy_sec_robust, policy_sec.dvresreform ~ policy_sec.cons_ca
                      , firth=TRUE)
  summary(sr1_test)
  #
  sec_lib_1c_test <- logistf(data=policy_sec_robust, policy_sec.dvlibreform ~ policy_sec.cons_ca
                             , firth=TRUE)
  summary(sec_lib_1c_test)
  
  texreg(list(sec_lib_1c_test, sr1_test),
       #   file="tableA5.2.pdf",
         custom.coef.names = c(
           "Intercept",
           "Strength Conservative Parties"
         ),
         reorder.coef = c(2,1),
         caption= "Secular World: with Italy and Ireland",
      #   label="tab:5",
         single.row=TRUE
      #   sideways=TRUE
  )  
  
####################################################################
  #descriptive statistics
  #Tables A6.1 and A6.2:
  library(stargazer)
  stargazer(subset(policy_sec, select=c(dvlibreform, dvresreform, cons_ca, v7210_pii, kat, loggdb, tsR2, i_all_perc)))
  
  stargazer(subset(policy_rel, select=c(dvlibreform, dvresreform, cons_ca, v7210_pii, kat, loggdb, tsR2, i_all_perc)))

######################################################################
  
  #Prepare the descriptive graphs to explore how strong conservative parties have been for policy reforms:
  
  library(beeswarm)
  
  policy_rel_res <- subset(policy_rel, dvresreform==1, select=cons_ca)
  policy_rel_res$reform <- 1
  names(policy_rel_res) <- c("cons", "reform")
  
  policy_rel_lib <- subset(policy_rel, dvlibreform==1, select=cons_ca)
  policy_rel_lib$reform <- 2
  names(policy_rel_lib) <- c("cons", "reform")
  
  policy_sec_res <- subset(policy_sec, dvresreform==1, select=cons_ca)
  policy_sec_res$reform <- 3
  names(policy_sec_res) <- c("cons", "reform")
  
  policy_sec_lib <- subset(policy_sec, dvlibreform==1, select=cons_ca)
  policy_sec_lib$reform <- 4
  names(policy_sec_lib) <- c("cons", "reform")
  #######
  # Figure A3: Strength of Conservative Parties for Reforms (Secular World)
 
  
  
  
  
   
  jpeg("figure_A3.jpg", width=6, height=4, units='in', res=600)

  par(mfrow=(c(1,2)))
  
  beeswarm(policy_sec_res$cons, data = policy_sec_res, 
           log = FALSE, pch = 16, col = "black",
           main = 'Restrictive Reforms', 
           ylab="share of conservative cabinet seats", 
           ylim=c(0, 1))

  beeswarm(policy_sec_lib$cons, data = policy_sec_lib, 
           log = FALSE, pch = 16, col = "black",
           main = 'Permissive Reforms', 
           ylab="share of conservative cabinet seats", 
           ylim=c(0, 1))

  dev.off()
  
  #####
  
  # Figure A4: Strength of Conservative Parties for Reforms (Religious World)
  
  jpeg("figure_A4.jpg", width=6, height=4, units='in', res=600)

    par(mfrow=(c(1,2)))
  
  
  beeswarm(policy_rel_res$cons, data = policy_rel_res, 
           log = FALSE, pch = 16, col = "black",
           main = 'Restrictive Reforms', 
           ylab="share of conservative cabinet seats", 
           ylim=c(0, 1))
  
  
  beeswarm(policy_rel_lib$cons, data = policy_rel_lib, 
           log = FALSE, pch = 16, col = "black",
           main = 'Permissive Reforms', 
           ylab="share of conservative cabinet seats", 
           ylim=c(0, 1))
  dev.off()
  
  ##############################
  
  #Figure A5: Religious World - liberal reform: by policy area
  
  policy_rel_lib_abo <- subset(policy_rel, (policy_rel_1==1 & libeventsabo==1), select=cons_ca)
  policy_rel_lib_euth <- subset(policy_rel, (policy_rel_2==1 & libeventseuth==1), select=cons_ca)
  policy_rel_lib_homo <- subset(policy_rel, (policy_rel_3==1 & libeventshomo==1), select=cons_ca)
  policy_rel_lib_pro <- subset(policy_rel, (policy_rel_4==1 & libeventspro==1), select=cons_ca)
  policy_rel_lib_porn <- subset(policy_rel, (policy_rel_5==1 & libeventsporn==1), select=cons_ca)
  policy_rel_lib_ssm <- subset(policy_rel, (policy_rel_6==1 & libeventsssm==1), select=cons_ca)

  jpeg("figure_A5.jpg", width=6, height=4, units='in', res=600)  
  par(mfrow=(c(1,6)))
  
  beeswarm(policy_rel_lib_abo$cons_ca, data = policy_rel_lib_abo, 
           log = FALSE, pch = 16, col = "black",
           main = 'Abortion', 
           ylab="share of conservative cabinet seats", 
           ylim=c(0, 1))
  
  beeswarm(policy_rel_lib_euth$cons_ca, data = policy_rel_lib_euth, 
           log = FALSE, pch = 16, col = "black",
           main = 'Euthanasia', 
           #     ylab="share of conservative cabinet seats", 
           ylim=c(0, 1))
  
  beeswarm(policy_rel_lib_homo$cons_ca, data = policy_rel_lib_homo, 
           log = FALSE, pch = 16, col = "black",
           main = 'Homosexuality', 
           #      ylab="share of conservative cabinet seats", 
           ylim=c(0, 1))
  
  beeswarm(policy_rel_lib_pro$cons_ca, data = policy_rel_lib_pro, 
           log = FALSE, pch = 16, col = "black",
           main = 'Prostitution', 
           #       ylab="share of conservative cabinet seats", 
           ylim=c(0, 1))
  
  beeswarm(policy_rel_lib_porn$cons_ca, data = policy_rel_lib_porn, 
           log = FALSE, pch = 16, col = "black",
           main = 'Pornography', 
           #        ylab="share of conservative cabinet seats", 
           ylim=c(0, 1))
  
  beeswarm(policy_rel_lib_ssm$cons_ca, data = policy_rel_lib_ssm, 
           log = FALSE, pch = 16, col = "black",
           main = 'SSM', 
           #         ylab="share of conservative cabinet seats", 
           ylim=c(0, 1))
  
  
  dev.off()
  
###################
  #Figure A6: Religious World - restrictive reform: by policy area

  policy_rel_res_abo <- subset(policy_rel, (policy_rel_1==1 & reseventsabo==1), select=cons_ca)
  policy_rel_res_euth <- subset(policy_rel, (policy_rel_2==1 & reseventseuth==1), select=cons_ca)
  policy_rel_res_homo <- subset(policy_rel, (policy_rel_3==1 & reseventshomo==1), select=cons_ca)
  policy_rel_res_pro <- subset(policy_rel, (policy_rel_4==1 & reseventspro==1), select=cons_ca)
  policy_rel_res_porn <- subset(policy_rel, (policy_rel_5==1 & reseventsporn==1), select=cons_ca)
  policy_rel_res_ssm <- subset(policy_rel, (policy_rel_6==1 & reseventsssm==1), select=cons_ca)
  
  
  jpeg("figure_A6.jpg", width=6, height=4, units='in', res=600)  
  par(mfrow=(c(1,6)))
  
  beeswarm(policy_rel_res_abo$cons_ca, data = policy_rel_res_abo, 
           log = FALSE, pch = 16, col = "black",
           main = 'Abortion', 
           ylab="share of conservative cabinet seats", 
           ylim=c(0, 1))
  
  
  beeswarm(policy_rel_res_euth$cons_ca, data = policy_rel_res_euth, 
           log = FALSE, pch = 16, col = "black",
           main = 'Euthanasia', 
           #     ylab="share of conservative cabinet seats", 
           ylim=c(0, 1))
  
  
  beeswarm(policy_rel_res_homo$cons_ca, data = policy_rel_res_homo, 
           log = FALSE, pch = 16, col = "black",
           main = 'Homosexuality', 
           #      ylab="share of conservative cabinet seats", 
           ylim=c(0, 1))
  
  beeswarm(policy_rel_res_pro$cons_ca, data = policy_rel_res_pro, 
           log = FALSE, pch = 16, col = "black",
           main = 'Prostitution', 
           #       ylab="share of conservative cabinet seats", 
           ylim=c(0, 1))
  
  beeswarm(policy_rel_res_porn$cons_ca, data = policy_rel_res_porn, 
           log = FALSE, pch = 16, col = "black",
           main = 'Pornography', 
           #        ylab="share of conservative cabinet seats", 
           ylim=c(0, 1))
  
  beeswarm(policy_rel_res_ssm$cons_ca, data = policy_rel_res_ssm, 
           log = FALSE, pch = 16, col = "black",
           main = 'SSM', 
           #         ylab="share of conservative cabinet seats", 
           ylim=c(0, 1))
  
  
  dev.off()
  
  ################################################
  #Figure A7: Secular World - liberal reform: by policy area
   
  policy_sec_lib_abo <- subset(policy_sec, (policy_sec_1==1 & libeventsabo==1), select=cons_ca)
  policy_sec_lib_euth <- subset(policy_sec, (policy_sec_2==1 & libeventseuth==1), select=cons_ca)
  policy_sec_lib_homo <- subset(policy_sec, (policy_sec_3==1 & libeventshomo==1), select=cons_ca)
  policy_sec_lib_pro <- subset(policy_sec, (policy_sec_4==1 & libeventspro==1), select=cons_ca)
  policy_sec_lib_porn <- subset(policy_sec, (policy_sec_5==1 & libeventsporn==1), select=cons_ca)
  policy_sec_lib_ssm <- subset(policy_sec, (policy_sec_6==1 & libeventsssm==1), select=cons_ca)

  jpeg("figure_A7.jpg", width=6, height=4, units='in', res=600)  
  par(mfrow=(c(1,6)))
  
  beeswarm(policy_sec_lib_abo$cons_ca, data = policy_sec_lib_abo, 
           log = FALSE, pch = 16, col = "black",
           main = 'Abortion', 
           ylab="share of conservative cabinet seats", 
           ylim=c(0, 1))
  
  beeswarm(policy_sec_lib_euth$cons_ca, data = policy_sec_lib_euth, 
           log = FALSE, pch = 16, col = "black",
           main = 'Euthanasia', 
           #     ylab="share of conservative cabinet seats", 
           ylim=c(0, 1))
  
  beeswarm(policy_sec_lib_homo$cons_ca, data = policy_sec_lib_homo, 
           log = FALSE, pch = 16, col = "black",
           main = 'Homosexuality', 
           #      ylab="share of conservative cabinet seats", 
           ylim=c(0, 1))
  
  beeswarm(policy_sec_lib_pro$cons_ca, data = policy_sec_lib_pro, 
           log = FALSE, pch = 16, col = "black",
           main = 'Prostitution', 
           #       ylab="share of conservative cabinet seats", 
           ylim=c(0, 1))
  
  beeswarm(policy_sec_lib_porn$cons_ca, data = policy_sec_lib_porn, 
           log = FALSE, pch = 16, col = "black",
           main = 'Pornography', 
           #        ylab="share of conservative cabinet seats", 
           ylim=c(0, 1))
  
  beeswarm(policy_sec_lib_ssm$cons_ca, data = policy_sec_lib_ssm, 
           log = FALSE, pch = 16, col = "black",
           main = 'SSM', 
           #         ylab="share of conservative cabinet seats", 
           ylim=c(0, 1))
  
  
  dev.off()
  
  ##########################################
  #Figure A8: Secular World - restrictive reform: by policy area
  
  policy_sec_res_abo <- subset(policy_sec, (policy_sec_1==1 & reseventsabo==1), select=cons_ca)
  policy_sec_res_euth <- subset(policy_sec, (policy_sec_2==1 & reseventseuth==1), select=cons_ca)
  policy_sec_res_homo <- subset(policy_sec, (policy_sec_3==1 & reseventshomo==1), select=cons_ca)
  policy_sec_res_pro <- subset(policy_sec, (policy_sec_4==1 & reseventspro==1), select=cons_ca)
  policy_sec_res_porn <- subset(policy_sec, (policy_sec_5==1 & reseventsporn==1), select=cons_ca)
  policy_sec_res_ssm <- subset(policy_sec, (policy_sec_6==1 & reseventsssm==1), select=cons_ca)

  
  jpeg("figure_A8.jpg", width=6, height=4, units='in', res=600)  
  par(mfrow=(c(1,6)))
  
  beeswarm(policy_sec_res_abo$cons_ca, data = policy_sec_res_abo, 
           log = FALSE, pch = 16, col = "black",
           main = 'Abortion', 
           ylab="share of conservative cabinet seats", 
           ylim=c(0, 1))
  #bxplot(cons_ca, data = policy_sec_res_abo, add = TRUE)
  
  
  beeswarm(policy_sec_res_euth$cons_ca, data = policy_sec_res_euth, 
           log = FALSE, pch = 16, col = "black",
           main = 'Euthanasia', 
           #     ylab="share of conservative cabinet seats", 
           ylim=c(0, 1))
  
  beeswarm(policy_sec_res_homo$cons_ca, data = policy_sec_res_homo, 
           log = FALSE, pch = 16, col = "black",
           main = 'Homosexuality', 
           #      ylab="share of conservative cabinet seats", 
           ylim=c(0, 1))
  
  
  beeswarm(policy_sec_res_pro$cons_ca, data = policy_sec_res_pro, 
           log = FALSE, pch = 16, col = "black",
           main = 'Prostitution', 
           #       ylab="share of conservative cabinet seats", 
           ylim=c(0, 1))

  beeswarm(policy_sec_res_porn$cons_ca, data = policy_sec_res_porn, 
           log = FALSE, pch = 16, col = "black",
           main = 'Pornography', 
           #        ylab="share of conservative cabinet seats", 
           ylim=c(0, 1))
  
  beeswarm(policy_sec_res_ssm$cons_ca, data = policy_sec_res_ssm, 
           log = FALSE, pch = 16, col = "black",
           main = 'SSM', 
           #         ylab="share of conservative cabinet seats", 
           ylim=c(0, 1))
  
  
  dev.off()
  
  ######
  ######
  #Sensitivity Analysis:
  
  # without abortion:
  policy_rel_noabo <- subset(policy_rel, policy_rel_1 ==0)
  policy_sec_noabo <- subset(policy_sec, policy_sec_1 ==0)
  
  sens_rellib1 <- logistf(data=policy_rel_noabo, dvlibreform ~ cons_ca  
                          , firth=TRUE, method="Penalized ML")
  
  summary(sens_rellib1)
  
  sens_relres1 <- logistf(data=policy_rel_noabo, dvresreform ~ cons_ca
                          , firth=TRUE, method="Penalized ML")
  summary(sens_relres1)
  
  
  sens_seclib1 <- logistf(data=policy_sec_noabo, dvresreform ~ cons_ca
                          , firth=TRUE, method="Penalized ML")
  summary(sens_seclib1)
  
  sens_secres1 <- logistf(data=policy_sec_noabo, dvlibreform ~ cons_ca 
                          , firth=TRUE, method="Penalized ML")
  summary(sens_secres1)
  
  
  
  ######
  #Without Abortion:
  library(forestplot)
  moral <- structure(list(
    mean = c(sens_seclib1$coefficients[2], sens_secres1$coefficients[2], sens_relres1$coefficients[2], sens_rellib1$coefficients[2]),
    lower = c(sens_seclib1$ci.lower[2], sens_secres1$ci.lower[2], sens_relres1$ci.lower[2], sens_rellib1$ci.lower[2]),
    upper = c(sens_seclib1$ci.upper[2], sens_secres1$ci.upper[2], sens_relres1$ci.upper[2], sens_rellib1$ci.upper[2])),
    .Names = c("mean", "lower", "upper"),
    row.names = c(NA,-4L),
    class ="data.frame")
  
  tabletext2 <- c("H1a: Secular World - Liberalization", "H1b: Secular World - Restrictive Reform", "H2: Religious World - Restrictive Reform", "H3: Religious World - Liberalization")
  
  jpeg("figure_A9_noabortion.jpg", width=6, height=4, units='in', res=600)
  
  forestplot(tabletext2, moral, new_page= FALSE,
             zero=0, lwd.zero=1,
             txt_gp=fpTxtGp(label = list(gpar(fontsize = 12)), ticks=gpar(fontsize=12, cex=1),
                            xlab=gpar(fontsize=12, cex=1)),
             col=fpColors(lines="black", box="black", zero="black"),
             boxsize=0.05, clip=c(-Inf,Inf), xticks=c(-5,-4,-3,-2,-1,0,1,2,3,4,5)
  )
  
  dev.off()
  
#####################
  ######No Euthanasia
  policy_rel_noeuth <- subset(policy_rel, policy_rel_2 ==0)
  policy_sec_noeuth <- subset(policy_sec, policy_sec_2 ==0)
  
  
  sens_rellib1_e <- logistf(data=policy_rel_noeuth, dvlibreform ~ cons_ca  
                            , firth=TRUE, method="Penalized ML")
  
  summary(sens_rellib1_e)
  
  sens_relres1_e <- logistf(data=policy_rel_noeuth, dvresreform ~ cons_ca
                          , firth=TRUE, method="Penalized ML")
  summary(sens_relres1_e)
  
  
  sens_seclib1_e <- logistf(data=policy_sec_noeuth, dvresreform ~ cons_ca
                          , firth=TRUE, method="Penalized ML")
  summary(sens_seclib1_e)
  
  sens_secres1_e <- logistf(data=policy_sec_noeuth, dvlibreform ~ cons_ca 
                          , firth=TRUE, method="Penalized ML")
  summary(sens_secres1_e)
  
  
  
  ######
  #Without Euthanasia:
  library(forestplot)
  moral <- structure(list(
    mean = c(sens_seclib1_e$coefficients[2], sens_secres1_e$coefficients[2], sens_relres1_e$coefficients[2], sens_rellib1_e$coefficients[2]),
    lower = c(sens_seclib1_e$ci.lower[2], sens_secres1_e$ci.lower[2], sens_relres1_e$ci.lower[2], sens_rellib1_e$ci.lower[2]),
    upper = c(sens_seclib1_e$ci.upper[2], sens_secres1_e$ci.upper[2], sens_relres1_e$ci.upper[2], sens_rellib1_e$ci.upper[2])),
    .Names = c("mean", "lower", "upper"),
    row.names = c(NA,-4L),
    class ="data.frame")
  
  tabletext2 <- c("H1a: Secular World - Liberalization", "H1b: Secular World - Restrictive Reform", "H2: Religious World - Restrictive Reform", "H3: Religious World - Liberalization")
  
  jpeg("figure_A10_noeuthansia.jpg", width=6, height=4, units='in', res=600)
  
  forestplot(tabletext2, moral, new_page= FALSE,
             zero=0, lwd.zero=1,
             txt_gp=fpTxtGp(label = list(gpar(fontsize = 12)), ticks=gpar(fontsize=12, cex=1),
                            xlab=gpar(fontsize=12, cex=1)),
             col=fpColors(lines="black", box="black", zero="black"),
             boxsize=0.05, clip=c(-Inf,Inf), xticks=c(-5,-4,-3,-2,-1,0,1,2,3,4,5)
  )
  
  dev.off()

  
  ######No Homosexuality
  policy_rel_nohomo <- subset(policy_rel, policy_rel_3 ==0)
  policy_sec_nohomo <- subset(policy_sec, policy_sec_3 ==0)
  
  
  sens_rellib1_h <- logistf(data=policy_rel_nohomo, dvlibreform ~ cons_ca  
                            , firth=TRUE, method="Penalized ML")
  
  summary(sens_rellib1_h)
  
  sens_relres1_h <- logistf(data=policy_rel_nohomo, dvresreform ~ cons_ca
                            , firth=TRUE, method="Penalized ML")
  summary(sens_relres1_h)
  
  
  sens_seclib1_h <- logistf(data=policy_sec_nohomo, dvresreform ~ cons_ca
                            , firth=TRUE, method="Penalized ML")
  summary(sens_seclib1_h)
  
  sens_secres1_h <- logistf(data=policy_sec_nohomo, dvlibreform ~ cons_ca 
                            , firth=TRUE, method="Penalized ML")
  summary(sens_secres1_h)
  
  
  
  ######
  #Without Homosexuality:
  library(forestplot)
  moral <- structure(list(
    mean = c(sens_seclib1_h$coefficients[2], sens_secres1_h$coefficients[2], sens_relres1_h$coefficients[2], sens_rellib1_h$coefficients[2]),
    lower = c(sens_seclib1_h$ci.lower[2], sens_secres1_h$ci.lower[2], sens_relres1_h$ci.lower[2], sens_rellib1_h$ci.lower[2]),
    upper = c(sens_seclib1_h$ci.upper[2], sens_secres1_h$ci.upper[2], sens_relres1_h$ci.upper[2], sens_rellib1_h$ci.upper[2])),
    .Names = c("mean", "lower", "upper"),
    row.names = c(NA,-4L),
    class ="data.frame")
  
  tabletext2 <- c("H1a: Secular World - Liberalization", "H1b: Secular World - Restrictive Reform", "H2: Religious World - Restrictive Reform", "H3: Religious World - Liberalization")
  
  jpeg("figure_A11_nohomo.jpg", width=6, height=4, units='in', res=600)
  
  forestplot(tabletext2, moral, new_page= FALSE,
             zero=0, lwd.zero=1,
             txt_gp=fpTxtGp(label = list(gpar(fontsize = 12)), ticks=gpar(fontsize=12, cex=1),
                            xlab=gpar(fontsize=12, cex=1)),
             col=fpColors(lines="black", box="black", zero="black"),
             boxsize=0.05, clip=c(-Inf,Inf), xticks=c(-5,-4,-3,-2,-1,0,1,2,3,4,5)
  )
  
  dev.off()

  ######No Pornography
  policy_rel_noporn <- subset(policy_rel, policy_rel_4 ==0)
  policy_sec_noporn <- subset(policy_sec, policy_sec_4 ==0)
  
  
  sens_rellib1_p <- logistf(data=policy_rel_noporn, dvlibreform ~ cons_ca  
                            , firth=TRUE, method="Penalized ML")
  
  summary(sens_rellib1_p)
  
  sens_relres1_p <- logistf(data=policy_rel_noporn, dvresreform ~ cons_ca
                            , firth=TRUE, method="Penalized ML")
  summary(sens_relres1_p)
  
  
  sens_seclib1_p <- logistf(data=policy_sec_noporn, dvresreform ~ cons_ca
                            , firth=TRUE, method="Penalized ML")
  summary(sens_seclib1_p)
  
  sens_secres1_p <- logistf(data=policy_sec_noporn, dvlibreform ~ cons_ca 
                            , firth=TRUE, method="Penalized ML")
  summary(sens_secres1_p)
  
  
  
  ######
  #Without Pornography:
  library(forestplot)
  moral <- structure(list(
    mean = c(sens_seclib1_p$coefficients[2], sens_secres1_p$coefficients[2], sens_relres1_p$coefficients[2], sens_rellib1_p$coefficients[2]),
    lower = c(sens_seclib1_p$ci.lower[2], sens_secres1_p$ci.lower[2], sens_relres1_p$ci.lower[2], sens_rellib1_p$ci.lower[2]),
    upper = c(sens_seclib1_p$ci.upper[2], sens_secres1_p$ci.upper[2], sens_relres1_p$ci.upper[2], sens_rellib1_p$ci.upper[2])),
    .Names = c("mean", "lower", "upper"),
    row.names = c(NA,-4L),
    class ="data.frame")
  
  tabletext2 <- c("H1a: Secular World - Liberalization", "H1b: Secular World - Restrictive Reform", "H2: Religious World - Restrictive Reform", "H3: Religious World - Liberalization")
  
  jpeg("figure_A12_noporn.jpg", width=6, height=4, units='in', res=600)
  
  forestplot(tabletext2, moral, new_page= FALSE,
             zero=0, lwd.zero=1,
             txt_gp=fpTxtGp(label = list(gpar(fontsize = 12)), ticks=gpar(fontsize=12, cex=1),
                            xlab=gpar(fontsize=12, cex=1)),
             col=fpColors(lines="black", box="black", zero="black"),
             boxsize=0.05, clip=c(-Inf,Inf), xticks=c(-5,-4,-3,-2,-1,0,1,2,3,4,5)
  )
  
  dev.off()
  
  ######No Prostitution
  policy_rel_noprost <- subset(policy_rel, policy_rel_5 ==0)
  policy_sec_noprost <- subset(policy_sec, policy_sec_5 ==0)
  
  
  sens_rellib1_pro <- logistf(data=policy_rel_noprost, dvlibreform ~ cons_ca  
                            , firth=TRUE, method="Penalized ML")
  
  summary(sens_rellib1_pro)
  
  sens_relres1_pro <- logistf(data=policy_rel_noprost, dvresreform ~ cons_ca
                            , firth=TRUE, method="Penalized ML")
  summary(sens_relres1_pro)
  
  
  sens_seclib1_pro <- logistf(data=policy_sec_noprost, dvresreform ~ cons_ca
                            , firth=TRUE, method="Penalized ML")
  summary(sens_seclib1_pro)
  
  sens_secres1_pro <- logistf(data=policy_sec_noprost, dvlibreform ~ cons_ca 
                            , firth=TRUE, method="Penalized ML")
  summary(sens_secres1_pro)
  
  
  
  ######
  #Without Prostitution:
  library(forestplot)
  moral <- structure(list(
    mean = c(sens_seclib1_pro$coefficients[2], sens_secres1_pro$coefficients[2], sens_relres1_pro$coefficients[2], sens_rellib1_pro$coefficients[2]),
    lower = c(sens_seclib1_pro$ci.lower[2], sens_secres1_pro$ci.lower[2], sens_relres1_pro$ci.lower[2], sens_rellib1_pro$ci.lower[2]),
    upper = c(sens_seclib1_pro$ci.upper[2], sens_secres1_pro$ci.upper[2], sens_relres1_pro$ci.upper[2], sens_rellib1_pro$ci.upper[2])),
    .Names = c("mean", "lower", "upper"),
    row.names = c(NA,-4L),
    class ="data.frame")
  
  tabletext2 <- c("H1a: Secular World - Liberalization", "H1b: Secular World - Restrictive Reform", "H2: Religious World - Restrictive Reform", "H3: Religious World - Liberalization")
  
  jpeg("figure_A13_noprost.jpg", width=6, height=4, units='in', res=600)
  
  forestplot(tabletext2, moral, new_page= FALSE,
             zero=0, lwd.zero=1,
             txt_gp=fpTxtGp(label = list(gpar(fontsize = 12)), ticks=gpar(fontsize=12, cex=1),
                            xlab=gpar(fontsize=12, cex=1)),
             col=fpColors(lines="black", box="black", zero="black"),
             boxsize=0.05, clip=c(-Inf,Inf), xticks=c(-5,-4,-3,-2,-1,0,1,2,3,4,5)
  )
  
  dev.off()

  ######No SSM
  policy_rel_nossm <- subset(policy_rel, policy_rel_6 ==0)
  policy_sec_nossm <- subset(policy_sec, policy_sec_6 ==0)
  
  
  sens_rellib1_ssm <- logistf(data=policy_rel_nossm, dvlibreform ~ cons_ca  
                              , firth=TRUE, method="Penalized ML")
  
  summary(sens_rellib1_ssm)
  
  sens_relres1_ssm <- logistf(data=policy_rel_nossm, dvresreform ~ cons_ca
                              , firth=TRUE, method="Penalized ML")
  summary(sens_relres1_ssm)
  
  
  sens_seclib1_ssm <- logistf(data=policy_sec_nossm, dvresreform ~ cons_ca
                              , firth=TRUE, method="Penalized ML")
  summary(sens_seclib1_ssm)
  
  sens_secres1_ssm <- logistf(data=policy_sec_nossm, dvlibreform ~ cons_ca 
                              , firth=TRUE, method="Penalized ML")
  summary(sens_secres1_ssm)
  
  
  
  ######
  #Without SSM:
  library(forestplot)
  moral <- structure(list(
    mean = c(sens_seclib1_ssm$coefficients[2], sens_secres1_ssm$coefficients[2], sens_relres1_ssm$coefficients[2], sens_rellib1_ssm$coefficients[2]),
    lower = c(sens_seclib1_ssm$ci.lower[2], sens_secres1_ssm$ci.lower[2], sens_relres1_ssm$ci.lower[2], sens_rellib1_ssm$ci.lower[2]),
    upper = c(sens_seclib1_ssm$ci.upper[2], sens_secres1_ssm$ci.upper[2], sens_relres1_ssm$ci.upper[2], sens_rellib1_ssm$ci.upper[2])),
    .Names = c("mean", "lower", "upper"),
    row.names = c(NA,-4L),
    class ="data.frame")
  
  tabletext2 <- c("H1a: Secular World - Liberalization", "H1b: Secular World - Restrictive Reform", "H2: Religious World - Restrictive Reform", "H3: Religious World - Liberalization")
  
  jpeg("figure_A14_nossm.jpg", width=6, height=4, units='in', res=600)
  
  forestplot(tabletext2, moral, new_page= FALSE,
             zero=0, lwd.zero=1,
             txt_gp=fpTxtGp(label = list(gpar(fontsize = 12)), ticks=gpar(fontsize=12, cex=1),
                            xlab=gpar(fontsize=12, cex=1)),
             col=fpColors(lines="black", box="black", zero="black"),
             boxsize=0.05, clip=c(-Inf,Inf), xticks=c(-5,-4,-3,-2,-1,0,1,2,3,4,5)
  )
  
  dev.off()
  
#Additional Regression table to assess the level of significance for the impact of conservatie parties on restrictive reforms in the religious world 
#with one additional level of significance:
  
  #table A7 (appendix)
  
  texreg(list(sens_rellib1_ssm),
         #   file="tableA5.2.pdf",
         custom.coef.names = c(
           "Intercept",
           "Strength Conservative Parties"
         ),
         reorder.coef = c(2,1),
         caption= "Religious World- Liberalization - without SSM",
         #   label="tab:5"
         stars = c(0.001, 0.01, 0.05, 0.1),
         single.row=TRUE,
            sideways=TRUE
  )  
  
  